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ABSTRACT 


This thesis analyzes the dynamic stability of positively buoyant submersibles. Six 
degree-of-freedom equations of motion are uscd to compule steady state behavior with 
motion restricted to the vertical plane. Steady state solutions are analyzed for various 
conditions of buoyancy including changes in (1) the amount of excess buoyancy, (2) the 
location of the center of buoyancy, (3) the location of the centcr of gravity, as well as (4) 
the deflection of bow and stern planes. The equations of motion are then linearizcd around 
these steady state solutions to predict dynamic responsc in the vertical plane. The stability 
of each solution is then determined by eigcn value analysis. The study then cxpands the 
analysis to include all six degrces of frecdom (i.e., include stability analysis in the 


horizontal plane). Finally, numerical integration mcthods are used to verify the results. 
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I. INTRODUCTION 


A. PROBLEM STATEMENT 

Controlling emergency ascent situations on submersible vehicles such as dive plane 
jam recovery is of concern to the U.S. Navy. In order to control such situations, one must 
first be able to predict the dynamic response of positively buoyant submersibles. 

Dynamic rcsponse equations of motion descibe the manuevcring characteristics of 
submersible vehicles for six degrees of ftrecdom. These equations assume constant 
coefficients for hydrodynamic forces and moments approximated by zero {frequency added 
mass and damping terms plus the quadratic terms for drag forces. The constant coetficients 
vary for cach vehicle and are dependent on such things as vehicle body shape, location and 
magnitude of vehicle weight, location and magnitude of vehicle buoyancy, position of bow 
and stern planes, position of rudder, vehicle speed, vehicle mass characteristics, vehicle 
hydrodynamic coefficients, propeller rpm and control surface inputs. This thesis uses the 
equations of motion and hydrodynamic coefficients for a submerged Mark IX Swimmer 
Delivery Vehicle (SDV) devcloped by Smith, Crane, and Summey [Reference 1:1 1-16,21- 
31] to forecast the dynamic stability of a submersible in a positively buoyant condition. 

This study begins by using the six equations of motion to compute the steady state 
behavior of a submersible vehicle with motion restricted to the vertical plane. The steady 
State solutions in the vertical plane are calculated for various conditions of buoyancy 
including changes in the amount of excess buoyancy, the location of the center of 
buoyancy, the location of the center of gravity, as well as the deflection of bow and stern 
planes. The SDV's cquations of motion are then linearizcd around these steady state 


solutions to predict dynamic response motion in the vertical plane for the various conditions 


of buoyancy. Scvcral solutions are computed and the stability of each solution is 
determined by eigen value analysis. The thesis then expands thc analysis to include all six 
degrecs of freedom (i.c. include stability analysis in the horizontal plane). Finally, 


numerical integration mcthods are used to verify the results. 


B. EQUATIONS OF MOTION 
The six degrec of frecdom equations of motion for the submersible vehicle shown in 
Appendix A were taken from Smith, Crane, and Summey [Reference 1:11-16]. 


“ 


Diftercntiation with respect to time is denotcd by a dot over the variable; e.g. u= a 


These equations are referenced to a right-hand orthogonal axis system fixed in the body 
(vehicle) as shown in Figure 1. Since these equations are in reference to a body fixed axis 
systcm, the Eulcr anglcs of pitch (8),roll (o), and yaw (") are used to specify orientation 
with respect to the inertial reference system. Thc rotation sequcnce for 9, 8 and ¥, and thc 
Euler angle rates for 0. Q and yp shown in Appendix B were taken from Smith, Crane, and 
Summcy [Reference 1:18-20]. Major variables and parameters as dcfined by Smith, 
Crane, and Summcy [Reference 1:7-10] arc given bclow : 

1. Dynamic Variables 

u,V,W - Lincar vclocity components of vehicle with respcct to origin of 


body axes system rclative to fluid. 


loner - Angular velocity components of vehicle with respect to body 
axes system relative to inertial refcrencc system. 


DS ERTE - Hydrodynamic force components along body axes. 


K,M,N - Hydrodynamic moment components along body axcs. 


2. Mass Distribution Parameters 


Mass of the tlooded vehicle, including the mass of the 
entrained fluid. 


Weight of the flooded vehicle, including the weight of the 
entrained fluid (= gm; where g is the accelcration of gravity). 


Displacement volume of the vehicle. 


Buoyancy force acting on the vehicle (=p g V ). This is 


independent of the inertial mass distribution of the submersible 
vehicle, including whether or not it 1s flooded. 


XGVqlg - Coordinates of the CG (center of gravity) in the body axis 
system (Figure 1). Thesc will depend on the mass 
distribution of the vehicle, including the mass of the entrained 
fluid. 


XpY¥pZp - Coordinates of the CB (center of buoyancy) in the body axis 
system (Figure 1). These are independent of the mass 


distribution system, but may vary with the addition or removal 
of external appendages. 


Moments of inertia about the body system axcs, including the 
entrained fluid. 


Products of inertia about the body system axes, including the 
entrained fluid. 


3. Remaining Parameters 


Mass density of fluid medium 


Reference length used to nondimensionalize the hydrodynamic 
coctficicnts. 


oS) 


b(x), A(x) - Width and height of vehicle in its xy and xz planes, 
respectively, at location x measured in the body axes system 
(Figure 1). These quantities are required in the integrals 
detining the erossflow forces and moments in the equations of 


motion, and are tabulated within the Steady State Computer 
Program (Appendix C). 


a ee - Sternplane, bowplane and rudder deflection angles in radians 
nose’ tail 


(Figure 1). 
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Il. SYSTEM SOLUTIONS IN THE VERTICAL PLANE 


A. GENERAL 

In the steady state condition, the submersible will have reached constant linear and 
angular velocities. Therefore, the body fixed linear accelerations (u , v, w) and the body 
fixed angular accelerations (p , q , r) will be zero. Similarly, the vehicle will have reached 
a constant angle of pitch (8) making its derivatives (8) equal to zero. Since this analysis is 
restricted to steady state solutions in the vertical plane, the angle of roll (9) and its 
derivative (9) will be zero (in Chapter IV, the case where the angle of roll (9) is 180 
degrees will be discussed during the numerical integration analysis). The angle of yaw (¥) 
and its derivative (?) will likewise be zero due to the vertical plane restriction. It should be 
noted that had this analysis not been restricted to the vertical plane, steady state yaw (“) 
would not necessarily be zero, thereby allowing the angular velocities (p, q, r) to be non- 
zero. However, since this analysis was restricted to the cases where 6, @ ; y , D, Grateien 
are all zero, the equations of motion for six degrees of freedom for the steady state 


condition can be reduced to: 
¢ Longitudinal (Surge) Equation of Motion: 
(W -B)sin@ =| Xyy v2 + Xww w2 + Xp. uvd; + uw( X95. ds + X wan Opl| 
+ w2 | Xasds Os Xabob dy” 4 Xdror 5,") + u2X prop 
¢ Lateral (Sway) Equation of Motion: 


- (W-B) cos 6 sing= | Yy uv + Yyw vw + Yo, u2 op 


Xnose ; 5 ane 
2 : | CDy h(x) (v) +Cpz b(x) (w) | Uer(x) dx 
tal 


: Normal (Heavc) Equation of Motion 


- (W-B) cos @ cos @ = Zw uw + Zyy v2 +02 (Zp. dg+Z pp dp) 


*nose F ae 
: CDy h(x) (v)& +Cpz b(x)(w)* | - a Os 
Xiail 
- Roll Equation of Motion: 

-(yG W - yp B)cos 8 cos 9 + (zGW - zBB) cos 6 sin 


= Lie LY oa vw! “= u2 Kprop 


- Pitch Equation of Motion: 
(xG W - xB B)cos 8 cos 9 + (2G W - zp B) sin @ = 


+ 
} 
! 
} 
| 
~ 


My uw + Myy v2 + u2 (M,. Os + May db] 
ese : 
+ CDy h(x) (v)@ +Cpz b(x) (w)* | —¥— x dx 
1 “Ue glx) 
tail 
: Yaw Equation of Motion: 
-(xG W-xBB)cos @ sing + (yG W - yp B) sin @ = 
| Ny uv + Nyw vw + Ng, u2 bo, + u2 Nprop 
Xnose - et 
- Cpy h(x) (v)j* + Cpz b(x) (w)* | —*— x dx 
_ 1 Dy ) Dz )( | Ues(x) 
tail 
B. CONDITIONS 
1. Defining Additional Terms 
a. Excess Buoyancy, OB 
Excess buoyancy is defined as 0B = B - W wherc B Is the submersible's 


total buoyancy and W is the submcrsible's total weight. 


b. Longitudunal Center of Buoyancy, xgp 


The longitudinal center of buoyancy 1s defined as xGp = XG - Xp where 
XG 1s the longitudinal center of gravity with respect to the body fixed axis and Xp Is the 
longitudinal center of buoyancy with respect to the body fixed axis. 
c. Vertical Center of Buoyancy, ZG RB 
The vertical center of buoyancy is defined as gp =Zq - Zp where Zc is 
the longitudinal center of gravity with respect to the body fixed axis and Zp is the 


longitudinal center of buoyancy with respect to the body fixed axis (zGp 1s assumed to be 


positive). 


2. Assumed Conditions 


a. Lateral Centers of Gravity, yo , and Buoyancy, yp 


The lateral center of gravity and the center of buoyancy are assumed to be 


on the same centerline plane (y, = yp = 0). 
b. Propeller Speed,n (revolutions per minute) 


The propeller speed 1s assumed to be zero (n = Q). 


and N 


c. Propeller Coefficients, K prop 


prop 
From Smith, Crane, and Summey [Reference 1:30], the propeller 


coefficients are zero (Kprop = Nprop = 9). 


C. REVISING THE EQUATIONS OF MOTION 


Using the term for vertical center of buoyancy (zGp), the expression (ZG W - zpB) 
may be written as (ZGRW - Zp OB). Similarly, using the term for longitudinal center of 
buoyancy (XGp). the expression (xW - XpB) may be written as (xGRW - XpOB). Also, 
the term u?Xprop may be written as: 


uXprop = u2Cpa{n?-1}-u2Cppq (Yeammanded 1 = CpoA?2n2 - Cpou2 


Uactual 


where A IS a constant. 


Since the shaft speed (n) Is zero, the expression may be further reduced to 
2X --C 2 Suhbstituti hese ex sions plus the te fi ery ere os 
u“Aprop pou*. substituting these expressions plus the term for excess buoyancy 
(OB) and the assumed conditions revises the equations of motion for the six degree of 


freedom system as follows: 


° Longitudinal (Surge) Equation of Motion: 
-OB sin @ =| Xyy v2 + Xwy w2 + Xo uvd; + uw( Xo d5 + Xap Op] 


| Y 7 . 2 ie ws x 2)\' 
a u2| OS Te aha Obs a oror el ) - Cpg ur 


¢ Lateral (Sway) Equation of Motion: 


OB cos 6 sing= Yy uv + Yyw vw + Yo. u2 dy 


Xnose p | 

- _Cxpy h(x) v- +Cpz b(x) w2 Cn dx 
; ‘ei / 

Xtail 

. Normal (Heave) Equation of Motion: 


OB cos 0 cos d= Ve uw + Zyy V2 + U2 (Z,.. dstZop, dp} 


>NOSE _ . | 
= | CDy h(x) Veo +Cpz b(x) w2 | Veil) 
Xtail 
: Roll Equation of Motion: 
(2G BW - zpdB) cos 6 sin o= Ky uv + Ky vw 


: Pitch Equation of Motion: 
(xGB W-xB 5B) cos 8 cos o + (2GB W - zB 6B) Sime 


L 


enoses a 
7 | CDy h(x) ve +CDz b(x) iS | Ues(x) x dx 
Xtail ct 


: Yaw Equation of Motion: 


(-xGR W + XB 6B) eos @ sin o= Ny Os Nvw MW oo No, u2 6, 


Xnose : 
: | rye ey Gack aye ee 
| ~Dy a Vex) 
Xtal 


These six equations only have five unknowns: u,v,w,9, and 9. Therefore, two of 
these equations must be dependent and additional conditions are required in order to make 


the number of equations equal the number of unknowns. 


D. ADDITIONAL CONDITIONS 

The next eondition to be applied to the vehicle is that the rudder will remain 
eenterlined, that is Or = O. Sinee the solutions of interest are those in whieh the vehicle 
remains within the vertical plane, it can be further speeified that the linear velocity in the 
transverse direction (v) equals zero. Recalling that the angle of roll (9) has been previously 
assumed to equal zero, the trigonometrie funetions of 9 can be identified as sin? equals zero 
and eos equals unity. Substituting these quantities back into the equations of motion, 
three of the six equations of motion (sway, roll, and yaw) yield trivial solutions. In 


addition, the eross-flow veloeity term (U.,) for the heave and pitch equations can be 


redueed to: 


cr _ 


Ucdx) = (v +xr?? + (w - xq??? 10.5 


=|w2)5 = 
since v, r, and q are zero. Furthermore, since Ch, is constant, it can be taken outside the 


integral. Therefore, the three remaining equations ean be wriiten as: 
¢ Longitudinal (Surge) Equation of Motion: 


. Krew 2+ uv ( Dy ea Sara dp] 
+ u2| Xgcos Os + AQpob Ob } -Cpo u 
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° Normal (Heave) Equation of Motion: 


“ Pitch Equation of Motion: 
(xGB W-xp OB} cos 0 + (ZGR W -zp OB} sin @ = 
. ‘NOSE 5 
My uw + u2(My- ds +Mop dp} Cpz b(x) x x dx 
tail 


E. VERTICAL PLANE EQUATIONS OF MOTION 


By defining the terms A,, as the | b(x) dx and x, as L| b(x) x dx, the three 


Xtatl ‘rail 


remaining equations defining motion in the vertical plane can be written as: 

¢ Longitudinal (Surge) Equation of Motion: 

-8B sin @ = Xww w2t uw( X95. ds +X wan dp) 

2| 2 r ‘| 2) 

l= A ee Og a Xshob Ob = Cpo us 
¢ Normal (Heave) Equation of Motion: 

OB cos 8 = Zy uw tu2(Z>_ dgtZa, dp} - CDz w wl Ay 
: Pitch Equation of Motion: 


(xGB W-xp 6B) cos 0 + (27GB W -ZzB 6B) sin 0 = 
My uw + u2(My_ ds +Map dop}t Cr, w iwi xaAy 


F. COMPUTER PROGRAM DEVELOPMENT 
Taking these three equations which describe motion in the vertical plane, solving the 


first two for sine 6 and cosine 0, respectively; and dividing all three through by u? yields: 


: Longitudinal (Surge) Equation of Motion: 
w )2, (Ww | 
sin@ __ 1 | Xww (“+ (%)( Xyas Os + Xwob Ob) | 


Ms es 2 | 
ur BB +( Xp555 85° + XppapOd-] Coo | 
: Normal (Heave) Equation of Motion: 
cos8_ 1. WwW - w |w w |w| A,.| 
u2 - 5B Ly u ar Wa Og+Zsp dp} Ci 2 | 


. Pitch Equation of Motion: 
ixGp W- xB 5B) © oe + (2GB W - zp OB) ae 
u 


W |W 
My _ te (My. Og PE vi dp}+ ORE one SOAP 


Now, defining the quantity a as w, and substituting w' back into the three 


equations: 
if w is positive: 


* Longitudinal (Surge) Equation of Motion: 
sin 1, Xww(w'P+ (w')( Xyog Os + Xyap Ob) 


; | Z 2 
: OB + | X5c6s Os + AQpob Ob - Coo 


° Normal (Heave) Equation of Motion: 
=e : = Dy (w') + (Zys ds+Zo Ob) -CDz(w'? Aw 


° Pitch Equation of Motion: 


—- Se 


(xGR W-xpoB (2GR W - zp OB 
My, (w') + (Ms. Os +Mop Op) + Celia a 


however, if w is negative: 


Longitudinal (Surge) Equation a Motion: 
| 
sin@__ 1  *wwlw P+ (w')( Xyas Os + Xwon ob) 


| 
| | 

2 Z Z 
BB + Xpc55 85° + Xppap Ob-] = Coo 
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Q Normal (Heave) Equation of Motion: 


Coe lw 1 0 2 One Caw Ad 
2 ge) [Za5 Os*+Zyp Ob) + CDz(w') 





: Pitch Equation of Motion: 
IXGB W-XB oB) se + (ZGB W -zB OB) wae = 


u 


Mw (w')+ (My. dg +Mo, db}- Cpz (wi? xaAy 


nee cos 8 


ee | sin cos 6 , , ee 
Substituting the equations for e and y2 into the pitch equation yields the 


— 


following expressions: 


if w is positive: 
(xGB W - xp dB) - 2s (ve) Zee ZeemOn = Cipy (we ae 


i— . 


oie + (wil X PR ORT tte) 
: (2GB W - zp OB) = es | a wos °S ; wodb b] 


Vig (w') + (M,. Oc + Mop dp} + Cpz (w'f SAAw 


and if w is negative: 


2 


(xGp W-xp OB) > Ws (wl 1Z on Zee On| 2 Cpcw Ay 


yy 1\? Whee : 
a4 Xww (wi i I Awads Os * Awab ob) 
Oe a | Xasas Os + ASbob ob} - Cho 


- (zop W - zp OB) 


My (w') + (Myo Os +Moy Op} - Caz (w'P xan 


Rearranging these two equations to get them into the form: A(w')? + B(w') + C = 0, the 


expressions become: 


if w is positive: 
CDz XAAw + (2GB W - 2B 6B) 3 Xww 


+(xGp W-xp 0B) 1 Cp, A, 
i) Mw + (2GB W -zp OB) 1 | Awos Os + Awad bb) 
6B 


; ( 
| -(xGB W- xp 0B) | Zy, 
| oB 


(w'f’ 





+ 


(zGR W - ZB 6B) a X55 0 Nee ot XSpdb aaa - can 


fs @ 
+ (Ms. Og + Mo, dp} -(xGB W - xB 6B) (28s do tZap, dp). 


| 


. 
"] 


' 


and if w is negative: 


| - C—z XAAw + (2GB W -zB 6B) ne Aww 
OB a 
-(xgB W-xp oB) 1 Cp, Ay 
My + (zGB W - zB OB) | Kee on: Xen Oe) | 
4 w| 
| -(xgp W -xp dB) z,, 
(2GB W - ZB 6B) = rane ne 9 XSbdb ane = C5 


+ (Myo Os +My, Ob) -(xGB W - xp OB) = (2s ds+Zop Ob) 


=i 


| 
| 
L 





+ 





These quadratic expressions were then solved using the equation: 
-B + | B*-4AC 
2A 








where: 
B =0B My +(zgp W-zp dB) ( X95. d5 + Xap Obl 


-(xGB W - xB OB) Zy 
C= (2gR W-zp OB) (Xp .9. O57 + Xppay dp~ - Coo] 
+ 0B(Ms. ds +My, dp)-(xGB W-xB dB) (Z. dg+Zo, dp} 


if w is positive: 


A = dB(Cp; xaAw) + (2GB W - 2p 8B) Xww +(xGBW-xBoB) Cpz As; 
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and if w is negative: 
A = - OB(Cpyz xaAy) + (26GB W - zB OB) Note (xo W-xp OB) Sinyecser 
The value of 8 was determined using the computed values of w' and the equation for 


tangent 0: 
sin 0 


u2 
tan 8 = —=: 
cos 0 


5 


uc 
ge0s Ge 
However, the value of ae varied depending on the value of w', which lead to two 
possible solutions: 


Equation for tan 0 if w' is Positive: 


: ee é Me 2 
-Xww OF -()( Xyvas Os + Xwab Ob) -( Xosos 2s” + Xopap Ob) + Coe 


tan Q = - - —-- —- 
Zyy (w') + + Ze dstZay, dp] -Cp, Afw']w'" 
Equation for tan 0 if w' is Negative: 
; nO 4 2 ‘ 2 y 2 
ems a SY (Xo oe Xwab Ob) -( Xosas OS memeN A aOD | + Cyo 


7 7! a 
Zn (W') + + (2% OotZoy dp| + Cp, Ay{w' jw 

In cither case, the value of u* was computed using the expression derived from the 
surge equation of motion: 


= jee, Oeste 


‘\, t m . r r 2 : a) 
- Xww (w'P -(w')[ Xe Os + Xap Ob} | X56 Os + AOhop Ob | + Coo 


This leads to two possible solutions for u (ic. u=+ \)u2). The value of w was 
computed using w = u (w’). Combining the two possible solutions of u with the two 
possible values for w' derived from the quadratic expression lead to four possible 
combinations of solutions for u and w. The computer program which calculates these four 
possible solutions is contained in Appendix C. It is an interactive program designed to 
allow the operator to select the amount of excess buoyancy as a percentage of vehicle 


weight (OB), the deflection of dive planes in degrees (0,), the ratio of bow planes to dive 


in 


I 


weight (OB), the deflection of dive planes in degrecs (0,), the ratio of bow planes to dive 
planes (0)/0,). the location of xGp and Xp from body fixed axis origin as a percentage of 


length, and the location of zap and zp, from body fixed axis origin in fect. 


G. STEADY STATE RESULTS 
Figures 2, 3, and 4 show typical steady state solutions for surge velocity, heave 


vclocity. and pitch angle as a function of dive plane angle. Thc two cases shown vary the 


location of the longitudinal center of buoyancy; for case (a): XG, = - | % of the vehicle 


B 
length (L), and for case (b): xgp=+ 1% L. The following parameters were the same for 
both cases: excess buoyancy , 0B = 2 % of the vehicle weight (W); deflection of bow 


planes, Op = 9; location of horizontal and vertical centers of buoyancy, x, = zp, = 0; and 


location of vertical center of gravity, Zp = 0.1 feet. 


All runs developed four solutions. For two of the solutions, the magnitude of the 
surge vclocititics were large while the magnitudes of the associated heave velocitics were 
rclatively small. This has been descibed by Booth [Ref 2: 297] as "predominantly forward 
motion". The other two solutions had small surge velocity magnitudes and large heave 
velocity magnitudes. Booth [Ref 3: 346] rcferred to this type of motion as "nearly vertical 
ascents". The positive or ncgative nature of the velocities are associated with the value of 
pitch angle. Positive heave velocitics are associated with pitch anglcs greater than 90 
degrees. That is thc submcrsible would be ascending in a belly up orientation. Although 
this steady state analysis computes four possible solutions, it gives no indication as to 
which of the solutions are stable (if any). Dynamic response and stability criteria will be 


discussed in the next chapter. 
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Figure 2. Steady State Vertical Plane Solutions for Surge Velocity (u) 
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Figure 3. Steady State Vertical Plane Solutions for Heave Velocity (w) 
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Figure 4. Steady State Vertical Plane Solutions for Pitch Angle (6) 
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HI. DYNAMIC STABILITY 


A. GENERAL 

The first portion of this chapter uses the six degrce of freedom equations of motion 
along with the Euler anglc rate equations for the derivatives of the angles of pitch and roll 
(0 and 6) to predict the dynamic stability when movemcnt is restricted to the vertical plane. 
These equations of motion are then linearized around the vertical plane steady state nominal 
points computed in chapter I. Eigen value analysis is then used to compute the stability of 
each solution. The second part of the chaptcr expands this analysis to include all six 
degrees of freedom and uses the same steady state nominal points to predict the dynamic 


stability when motion is not restricted to the vertical plane. 


B. RESTRICTING MOTION TO THE VERTICAL PLANE 

Since motion is restricted to the vertical plane, thc body-fixed transverse veloctiy (v) 
and its dcrivative (v) are zero. The angles of roll (@) and yaw (¥ ), and their derivatives (9 
and Y) are also zero. We will continue to assume that the lateral centcr of gravity and the 


lateral center of buoyancy arc on the same centerline plane (yg = yp = Q), and the rudder ts 


centerlined (Or =Q). 


C. LINEARIZED VERTICAL PLANE EQUATIONS OF MOTION 
Substituting these values into the original equations of motion (Appendix A), yields 
trivial results for three of the six equations: Lateral (sway), Roll, and Yaw. The remaining 


equations reducc to the following form: 
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¢ Longitudinal (Surge) Equation of Motion: 
Beet 4 ae 2 2, 9 
Pollle= De Use 2G G4 =| -Cpv - Naencs +Xapopob ag 
= i Xy9.08 + XyapOb: UW SE X05 > Xg5pOD. uqgt Pre w- 


+ Xwq> Mm wq+§Xoq + mxg! q2 - (W - B) sin B 


° Normal (Heave) Equation of Motion: 
m-Zy W-im xgt Z, gq =(Zp,08 + Zp, 0b! u? 


+ Z, uw + om ug+ mz q? 


ynose 


- Cc 7 b(x)(w-x 2A 2) dx + (W - B) cos 8 
bz b(x)(w-xq)- Tee) ( ) 


a»tail 


. Pitch Equation of Motion: 


imzg| u-.M,, +m xg] w + | Ty - Mg'q =,M,,s + Mg,0b; u2+fMyi uw 


Xnose 


2, (w - xq) 
evi. = mx ; Ud- MZG Wd - C 7 DOOvexG)=(- = = x dx 
Mg G uq- MZG wq Cz bOx)(w-xq) Uats) 


stail 


< (xGWw e xpB) cos 6 - (Z¢,W = ZpB) sin 6 


These three equations which deseribe motion restricted to the vertical plane are functions of 
four variables (u,w,q,@) plus their deriatives (u,w,q,6). The equations of motion and the 


Euler angle rate equation for 8 were linearized using the following generalized procedure: 


af; . afi | of, ] Of 
ete a wy qte—) 


B11(i,1)0 + BI(i,2)w + B11(i,3)g + B11(i,4)0 = um {OW Led’ 


0 


| * 
| 
A 


d9 @, 


where the B11(i)'s are the constants associated with the derivatives of the variables, the 


functions f; represents the right hand side of the nonlinear equations, and 1 = 1 to 4 


Za 


identifies the three equations of motion (surge, heave, and yaw) plus the Euler angle rate 


equation for 8. The partial derivatives were computed as follows: 


. Partial Derivatives of the Longitudinal (Surge) EOM: 


dt 


au = “lowe =e Delete Xe eabain aF (X 0S a X yapdbhwo = A11(1,1) 


ty _ Ji NG + (x OS + X opdb}ug = A11(€1,2) 


dw wos 


- = (Xwq - m)wot [X goss + X gondb}up = Adi (le) 


dh = (WW - B) cos Qo = ALI(1.4) 
a0 


° Partial Derivatives of the Normal (Heave) EOM: 


- = ZW 0) +2 (Z,.08 + Zn, Ob)ug = Al1(2,1) 


df> 
o2 _ ype A wee 
Ow 


as o(Zq + m) + 2Cp,Ayxalwol = A11(2,3) 


OP = (W-B) sin 09 = Al1(2,4) 
a0 


Note on differentiation procedure: The cross-flow velocity term (U,) was 


2h eae eee 


reduced to Uedx) = (v +xry° + (w - xq)! =|w -q|. Allowing 


Z2 


the integral term in the heave equation to be defined as I. where: 


le= co b(x)(w-xq)w - xq| dx 
xtail 





al, S¢ | xnoase 
— = 2Cp. b(; -xq)s = dx = 2Cn. b(x) lw - xal 
aaa »f | (x)(w-xq)sign (w - xq) dx >| (x) |w - xqi dx 
dl al d(w - a 
aa {| aw a Fe 2 = 2Colvil b(x) dx = 2Cp, Wo! Aw 
xtail 


al al, \{ a(w- ~ 
“ | ene oa D| = - 2Cp hwo | i b(x) dx = - 2Cp, |wol xa Aw 


Partial Derivatives of the Pitch EOM: 


3 = 2(M,.ds + M,,ob}ug + My wo = A11(3.1) 

u 

dt , “~ 

i My, ug +2ugCp,AyX Awol = A11(3.2) 

aq = (M, - mxg)ug -(MZG)Wo - Zia I Alwo| = Al 33) 


$f5 _ (xgW - xpB) sin®o - (2qW -zpB) cos 09 = ALI(3.4) 
a0 


Note on differentiation procedure: Again the cross-flow velocity term (U_,) was 


reduced to |w - q|. Allowing the integral term in the pitch equation to be defined 


as L,, where: 


In = co b(x)(w-xq)w - xq| x dx 


aq nos 
a(w 7) = xo b(x) lw = xq x dx 


ail 


xnose 


= 2Cp1Wo) | b(x) x dx = 2Cp, [wo] Xa Aw 


ae 


OW 





Ow 


XNOSE 


oF | al; | ay i v = cova x2 b(x) d= = IR iz lwo IA 





aq \a(w-q) aq 


. Partial Derivatives of the Euler Angle Rate Equation for Q: 


at, 


oa) ane 
du (4,1) 
or cle 

Ow 

at 


a= = iis 
. (4,3) 


20 Alia) 


at 
ele) 


The constants corresponding to the derivatives of the four variables (i.c. the left hand 


side of the four equations) are as follows: 
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Constants Associated with Derivatives for the Longitudinal 
(Surge) EOM: 


Mee Po = Ce) 


gq = mz, = BII(1,3) 


Blic.2) = Bi1(1,4) =0 


Constants Associated with Derivatives for the Normal 


(Heave) EOM: 


w = m-Zy = BI1(2,2 


q =-(mxg + Z,)= B11(2,3) 


B11(2.1) = B11(2,4) = 0 


Constants Associated with Derivatives for the Pitch EOM: 


y= -aize = Bilt) 


w => -(My +m xg) = BIIG,2 


q = Iy-Mj=Bl1G.3 


B11(3,4) = 0 


¢ Constants Associated with Derivatives for the Euler Angle Rate 


Equation for 6: 
q Sena 


B11(4,1) = B11(4,2) = B11(4,3) =0 
These expressions can be arranged tn a matrix format to form the linearized equations 
of motion in the vertical plane about the nominal steady state points. The matrix format 1s 


as follows: 


Bit« Xl=]All- Xa 


where: 


All(i,1)  =Al11(1,2) = A11(4,3) ~— A114) 

Pe ALI Wie ANI 2), eA) ee ne 
AVG) eA G2) el) ener) 
All(4,1) = A11(4,2) — A11(4,3) — A11(4,4) 
Bll(1,1) Bli(1,2) Bli(d,3)  Blid,4) 7 
es Bli@,1) “\BIi22)— (BIG 3 yee ners, 

| =BliG,!1)  B11(3,2) B11,3) B11@,4) 
B11(4,1)  Bll(4,2) B11(4,3) ~=B11(4,4) | 
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D. VERTICAL PLANE COMPUTER PROGRAM DEVELOPMENT 

The matrix format of the linearized equations of motion was used in the eomputer 
program shown in Appendix D to predict the dynamie stability of the vehicle with motion 
restricted to the vertical plane. The program is interactive in that it allows the operator to 
select which of the four data files from the steady state analysis (Chapter II) will be used to 
define the nominal points for the linearization process. An eigen system subroutine was 
used to find cigen values and eigen vectors. The program's output was called the degree of 
stability and only shows the largest real part of all eigen values. The stability criteria is 


such that the degree of stability must be negative in order for the solution to be stable. 


E. VERTICAL PLANE DYNAMIC RESULTS 

Of the four possible stcady state solutions computed in Chapter I, only one solution 
from each case yielded stable characteristics. There were some cases in which none of the 
solutions were stable for certain ranges of parametcrs. The general trend of the linearized 
dynamic results are fairly well demonstrated by the two cases discussed in Chapter IL. 


Recalling the parameters of these cases: 0B = 2%, ratio of bow to stern planes (6,/0.) = 0. 


Zcp = 0.1 fect, and x, =z, = 0. The first case placed the longitudinal center of gravity aft 
of the longitudinal center of buoyaney (Xg¢p = - | %), and the seeond ease placed the 


longitudinal center of gravity forward (xgp = + 1 %). Once again dive plane deflection 
angle (0.) was varied from - 20 to + 20 degrees for both cases. Figure 5 shows 


longitudinal velocity (u) as a function of dive planc angle (0,). Case One (X¢p = - | %) 


pad 


showed predominantly forward motion, while case two (XGp = - 1 %) yielded nearly 


vertical ascents. Figure 6 shows vertical motion (w) as as a function of dive plane angle 
(0.). The results concur with Figure 5, case one shows very little vertical motion while 
case two demonstrates a larger value. It is interesting to note that vertical motion (w) for 
case one Is positive for dive plane angles between - 20 and -4 degrees. The case one values 
of 8 shown in figure 7 for dive plane angles between - 20 and -4 degrees concur with this 
observation. The values of 6 greater than 90 degrees indicate the submersible is ascending 
in the belly up position. The stability in the vertical plane is shown in figure 8. Degree of 


stability (€) is shown as a function of dive plane angle (0.). This figure shows both cases 


to be stable in the vertical plane. 


u 
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Figure 5. Stable Vertical Plane Solutions for Surge Velocity 
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Figure 7. Stable Vertical Plane Solutions for Pitch Angle (6) 
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Figure 8 Degree of Stability in the Vertical Plane 


F. LINEARIZATION OF FULL EQUATIONS OF MOTION 
Referring back to the complete equations of motions shown in Appendix A, these six 
equations which describe motion for the six degree of freedom system are functions of 
eight variables (u,v.w,p.q.1,9,8) plus their deriatives Gee i). The sim 
equations of motion (Appendix A) along with the Euler angle rates (Appendix B) were 
linearized as folows: 
bjlu + bj2v + bj3w + bj4p + bj5q + bj6r + bj7o + bj86 = 


(ay z fe 7 re fe re a 
do. ldo: | ; d =| dg: 
OSI) yy | SBI ees A! eu eo r+ jnel] o 4] oe S) 


f du JUo L av IVo LOW jw L dp JPo L oq JQo or Ire [ ad Ido ee 410. 
where the bj's are the constants associated with the derivatives of the variables, the 


functions g; represents the right hand side of the nonlinear equations, and j = 1 to 8 


identifies the six equations of motion plus the Euler angle rate equations for 9 and 0. 


respectively. The partial derivatives of these equations were computed as follows: 


° Partial Derivatives of Longitudinal (Surge) EOM: 


OY : Gee | ( ; oe ; 
Oe wohl eerie Os GX nee obtuse pany Saul 


du 

dg ; : 

_ = 2XwwW aE (Xn OS + ais dbuy - 4913 

Og | ; 
a = -MWo + XwaW 0 nk xe OS + Xb db}ug - 415 


Bie -(W -B)cos 8 =al8 


a0 

VO r dy = 
°8l == al? logealg “EL s0=al6 SUE 
dv dp dr dg 


. Partial Derivatives of Lateral (Sway) EOM: 


dg) , 
vs =Yyug +} vwW() - Cp, Ay Wo = dle 
av 
I> 
“P= =mwo t+ Yo ug + YwpWo = a24 
op 
ago ie Lea 
ny = — ING Tee Ut Ywr yon CprAwXAWol= a26 
dy 
mam = (W - B) cos0 = a2/ 
dp 


: , d 
dg dg 7 ¥- o2 

£2 - Q=a2l C82 — (= a23 oB2 — Q) = 925 = = Q= a28 
du OW aq a0 


Note on differentiation procedure: The cross-flow velocity term (U_,,) Is given 


0.5 


by: Ucdx) = i(v +xr)? + (w - xq), . The integral term in the sway equation 


was defined as I.. where: 


Xnose _ : 
l= Coy h(x) (v+xr)@ +Cp, b(x) (w-xq)* | ¥ =X" dx 





























. Uer{x) 
Xtail 
Xnose 
= (I) Uatx)| dx 
; cf{x) 
Xtal 
[ 0 
a Pr , | Uer -(v + xr) i 
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Partial Derivatives of Normal (Heave) EOM: 
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Note on differentiation procedure: The integral term in the heave equation 


was defined as ie where: 
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Partial Derivatives of Pitch EOM: 
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Note on differentiation procedure: The integral term in the pitch equation 


was defined as I., where: 
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Note on differentiation procedure: The integral term in the yaw equation 


was defined as Re where: 
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° Partial Derivatives of Euler Angle Rate for 9: 
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The constants corresponding to the derivatives of the cight variables (1c. the left hand 


side of the cight equations) are as follows: 


Loe | 
~) 


Constants Associated with Derivatives for the Longitudinal 
(Surge) EOM: 


u = (m-X,)=b11 Qo — (mz) = bs 


bl2 = bi3— bi4 = blo — bly ble) 


Constants Associated with Derivatives for the Lateral (Sway) EOM: 
V ae (m = xo} = bz. p => (-mZG = Xe) = b24 
r => (mxg- Y;)=b26 


Constants Associated with Derivatives for the Normal (Heave) 
EOM: 


Woo (mia Zi) — bss q = (-mxg- Z)=b35 
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Constants Associated with Derivatives for the Roll EOM: 
v => (-mzg - Ky) = b42 r = (-K;)=b46 
p => (;,- K5) = b44 
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Constants Associated with Derivatives for the Pitch EOM: 
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Constants Associated with Derivatives for the Yaw EOM: 
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¢ Constants Associated with Derivatives for the Euler Angle Rate 
Equation for 9: 
gq => 12=b77 


Pak Di Dis le = iO] bso — OU 
¢ Constants Associated with the Derivatives for the Euler Angle 
Rate Equation for 0 : 
6 = al oisro 
Pole= DO2 =002 — Do) — Dofe— BS) = boo = Daz — 0 
These expressions can be arranged in a matrix format to define the dynamic equations 
of motion for the six degree of freedom system linearized about the nominal steady state 


points. The matrix format is B x X =A x X, where: 
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b18 
b28 
b38 
b48 
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b68 
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— © = © se Ce > 


S 


Ba eS =) = = 


These matrices may be reordered such that they will be of the form: 


Bll O 
O Ber? 


X1 |_fAll 0 | 


This ts accomplished by rewriting the X matrix such that X1 is the same matix used 


in the vertical plane analysis: 
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The A matrix is restructured into a matrix containing four 4 x 4 matrices with 


the Al2 and A21 matrices containing only the zero element. 
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Similarly, B is restructured into a matrix containing four 4 x 4 matrices with 


the B12 and B21 matrices containing only the zero element. 
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As discussed, the X1 matrix established to descibe the linearized dynamic stabilty in the 


vertical plane is the same as the X1 matrix within X. 


In addition, the All and B11 


matrices from the vertical plane analysis are also identical to those in X . That is: 
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The A22, B22, and X2 matrices represent the additional equations necessary to describe 


the linearized motion for all six degrees of freedom; henceforth referred to as the horizontal 


plane contributions. The eigen funetion for the six degree of freedom model is computed 


by taking the determinant as follows: 


5 


Alt -AB11 0 | : 
det | =0 =>(det!Al] -ABI1, }(det A22-2B22. )=0 
| 0 A22 - B22 | 
Sinee the eigen function will be the produet of these two determinants, the 
resulting eigen values will merely be the union of the vertical plane eigen values and the 
horizontal plane cigen values. The signifieanee of reducing the eigen value ealeulation from 
an 8 by 8 matrix problem to two 4 by 4 matrix problems is not in the computation time 


saved. But rather in faet that now the horizontal and vertical stabilities have been separated 


and identified. 


G. COMPUTER PROGRAM DEVELOPMENT 

The matrix format of the linearized dynamie response equations associated with 
motion in the horizontal plane was added to the eomputer program developed previously 
(Appendix D). Onee again, the program is interaetive in that it allows the operator to select 
whieh of the four data files from the steady state analysis (Chapter II) will be used to 
define the nominal points for the linearization proeess. An eigen system subroutine was 
used to find eigen values and eigen veetors. Two outputs were added to the program. 
First, the degree of stability in the horizontal plane, and next the degree of stability of both 
planes (i.e., the union of the vertical and horizontal degrees of stability). Reminder: 


degree of stability must be negative in order for the solution to be stable. 


H. DYNAMIC STABILITY SOLUTIONS 


Continuing with the same eases from part E, the horizontal stability for the two eases 


are shown in figure 9. Case two (xGp = + 1 %) Is stable in the horizontal plane for all 


values of dive planc angle (0). Whcreas, case one (Xgp = - 1 %) is unstable in the 
horizontal plane for dive plane angle (0,) between -20 and - 9 degrees. From figure 7, this 
corresponds to values of 0 grcater than 140 degrecs. This indicatcs that the vehicle is 
stable (even in the horizontal plane) for values of 8 greater than 90 degrees. A submersible 
with an angle of pitch greater than 90 degrees will have a negative metacentric height and 
will therefore be statically unstable. However, the results shown in figures 7 and 9 indicate 
that the vehicle will remain dynamically stable is such a condition. This ‘inverted 
pendulum’ type stability will be further investigated during the numerical integration 
analysis (Chapter IV) to see if hydrodynamic and drag forces on the vehicle can actually 
cause this to occur. Figure 10 shows the combined stabilities for the horizontal and vertical 
planes. It should be noted that in general the horizontal plane dictatcd stabilty for the cases 


considered. 
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Figure 9. Degree of Stability in the Horizontal Plane 
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Figure 10. Degree of Stability in Both (Horizontal and Vertical) Planes 
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Because the real part of the computed eigen values must be negative in order for a 
solution to be stable, it is easy to identify unstable solutions. However, the measure of 
stability for a stable solution is harder to quantify. The magnitudes shown thus far for 
degree of stability seem very small. Does this mean that the solutions are not very stable? 
In an attempt to answer this question, the values for degree of stability for a neutrally 
buoyant submersible (6B = 0 or W = B) were computed. Figure 11 shows the degree of 
stability in (a) the horizontal plane and (b) the vertical plane as a function of surge velocity. 
This figure shows that the magnitude of the values for degree of stability for this neutrally 


buoyant case are indeed of the same order of magnitude as the positively buoyant cases 


discussed earlier. 
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Figure 11. Degree of Stability as a Function of Surge Velocity (u) for a 


Neutrally Buoyant Submersible 
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IV. SIMULATIONS AND DISCUSSION OF RESULTS 


A. DYNAMIC STABILITY ANALYSIS RESULTS 

The dynamic stability analysis ineluded in this Chapter considers stability for both 
planes (1.c. the analysis no longer breaks stability down by horizontal or vertical plane). 
As mentioned in Chapter II, horizontal plane stability generally dictates the stability of the 
vehicle. 

1. Variations in Longitudinal Center of Gravity (Xe) 
Figures 12 through 15 show how ehanging the longitudinal eenter of gravity 


(X¢p) effeets the dynamic response of the submersible. For these eases, the amount of 
excess buoyancy (0B) ts two percent of wieght (W), the bow plane deflection angle (0,) Is 
zero, the vertical center of gravity (Z,,) 1s 0.1 feet, and the longitudinal and vertical eenters 
of buoyancy (x, and z,) are zero. The longitudinal center of gravity (xp) Is varied from 
- 2 to + 2 pereent of vehicle length. When the longitudinal center of gravity (x¢p) Is greater 
than zero, there are stable solutions for the full range of dive plane angles. However, this 
is not true when the longitudinal eenter of gravity (xGp) Is less than zero. The range of 


stable solutions is restricted when Xop=- O>eaind = 10: 
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Figure 12. 
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(b) xcp = 9, - 0.5, - 1.0, - 1.5 and - 2.0 % 


Stable Surge Velocity (u) Solutions for Variations in xcp 
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(a) xcp = 0, 0.5, 1.0, 1.8 and 2.0 % 
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(b) xcp = 9, - 0.8, - 1.0, - 1.8 and - 2.0 % 


Figure 13. Stable Heave Velocity (w) Solutions for Variations in x¢op 
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(b) xcp = 0, - 0.5, - 1.0, - 1.5 and - 2.0 % 


Stable Angle of Pitch (6) Solutions for Variations in xcp 
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Figure 15. Degree of Stablity for Variations in xox 
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2. Variations in the Amount of Excess Buoyancy (6B) 
Figures 16 through 19 show how changing the amount of excess buoyancy 


(OB) effects the dynamic response of the submersible. For these cases, the bow plane 


deflection angle (5,) is zero, the longitudinal center of gravity (XGp) Is - 0.5 percent of 
vehicle length, the vertical center of gravity (Z,,) 1s 0.1 feet, and the longitudinal and 
vertical centers of buoyancy (x, and z,) are zero. The amount of excess buoyancy (SB) is 
varied from 0.5 to 2.9 percent of weight (W). The only case where there are stable 


solutions for the full range of dive plane angles is when buoyancy (6B) is 0.5. 
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(b) OB = 2.5, 2.6, 2.7, 2.8, 29 & 


Figure 16. Stable Surge Velocity (u) Solutions for Variations in OB 
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(b) 6B = 2.5, 2.6, 2.7, 2.8, 2.9 % 


Figure 17. Stable Heave Velocity (w) Solutions for Variations in OB 
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Figure 18. 
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(b) 5B = 2.5, 2.6, 2.7, 2.8, 2.9 


Stable Angle of Pitch (6) Solutions for Variations in OB 
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(a) 6B = 0.5, 1.0, 1.5 and 2.0 % 
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(b) 6B = 2.5, 2.6, 2.7, 2.8, 2.9 % 


Figure 19. Degree of Stablity for Variations in 6B 
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3. Variations in Vertical Center of Gravity (Z-p) 


Figures 20 through 23 show how changing the vertical center of gravity (Zon) 
effects the dvnamic response of the submersible. For these cases. the amount of excess 
buoyancy (OB) ts two percent of weight (W), the bow plane deflection angle (0,,) is zero, 
the Jongitudinal] center of gravity (Xop) is - 0.5 percent of vehicle length, and the 
longitudinal and vertical centers of buoyancy (Xp and z,) are zero. The vertical center of 
gravity (Zp) is varied from 0.05 to 0.30 feet. The general trend shown in these figures is 


that the larger values for vertical center of gravity (Zp) have stable solutions over a larger 


range of dive plane angles. 
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Figure 20. Stable Surge Velocity (u) Solutions for Variations in Zcp 
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Figure 21. Stable Heave Velocity (w) Solutions for Variations in Zon 
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Figure 22. Stable Angle of Pitch (6) Solutions for Variations in z., 
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Figure 23. Degree of Stability for Variations in ZcR 


4. Variations in the Longitudinal Center of Buovancy (x;) 

Figures 24 through 27 show how changing the longitudinal center of buoyancy 
(Xp) effects the dynamic response of the submersible. For these cases, the amount of 
excess buovancy (OB) is two percent of weight (W). the bow plane deflection angle (0,) Is 
zero. the longitudinal center of gravity (Xp) is - 0.5 percent of vehicle length. the vertical 
center of gravity (Zon) is 0.1 feet. and the vertical center of buoyancy (Zp) Is zero. The 
longitudinal center of buovancy (Xp) 1s varied from - 9 to + 2 percent of vehicle length. 
The general trend shown in these figures is that the positive values for longitudinal center 
of buoyancy (Xp) tend to have stable solutions for positive dive plane angles. On the other 


hand, the negative values for longitudinal center of buoyancy (Xp) tend to have stable 


solutions for negative dive plane angles. 
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Figure 24. Stable Surge Velocity (u) Solutions for Variations in x), 
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Figure 25. Stable Heave Velocity (w) Solutions for Variations in Xp 
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Figure 26. Stable Angle of Pitch (6) Solutions for Variations in Xp 
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Figure 27. Degree of Stability for Variations in Xp 
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5. Variations in Bow Planes Deflection Angle (0,) 

Figures 28 through 31 show how a non-zero bow planes deflection angle (0.) 
effects the dynamic response of the submersible. For these cases, the amount of excess 
buovancy (OB) is two percent of wieght (W), the longitudinal center of gravity (Xp) Is 
- 0.5 percent of vehicle length, the vertical center of gravity (Z,,) is 0.1 feet, and the 
longitudinal and vertical centers of buoyancy (x, and z,) are zero. The bow planes are 


given a detlection of - 20 degrees. The significance of the results shown in these figures is 


that for certain dive plane angles (0. = -5 to -12 degrees) there are two stable solutions. 
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Figure 28. Stable Surge Velocity (u) Solutions for a Non-zero Bow Plane 
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Figure 29. Stable Heave Velocity (w) Solutions for a Non-zero Bow Plane 
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Figure 30. Stable Pitch Angle (6) Solutions for a Non-zero Bow Plane 
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Figure 31. Degree of Stability for a Non-zero Bow Plane Angle 


(0, = - 20 degrees) 


B. SIMULATIONS USING NUMERICAL INTEGRATION METHODS 
The linearized dynamic response results were verified by simulations using numerical 

integration of the full six degrees of freedom equations of motion for the swimmer delivery 

vehicle (SDV). Figure 32 shows a plot of angle of pitch (6) versus time for the center of 


gravity forward of the center of buoyancy case (Xgp = + 1) with a dive plane angle (0,) of 


- 15 degrees. The dotted line shows the linearized results from figure 7, the solid line 





shows the numerical integration results. The steady state results of the numerical 


integration method match the linearized dynamic results exactly. 
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Figure 32. Numerical Integration Solution for Angle of Pitch (6) when 


Center of Gravity is Forward (x, = + 1%) 
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Figure 53 shows a plot of angle of pitch (8) versus time for the center of gravity aft 
of the center of buoyancy case (X¢p = - 1) with a dive plane angle (0,) of - 15 degrees. 
Again the dotted line shows the linearized results from figure 7, the solid line shows the 
numerical integration results. And once again, the steady state results of the numerical 
integration method match the linearized dynamic results exactly. However, this linearized 
dynamic result was for the vertical plane only, the horizontal plane stability analysis 
indicated that this would be an unstable solution (figure 9). The reason for this 
disagreement in the results is investigated by adding an initial angle of roll to the numerical 
integration analysis. Adding a small angle of initial roll (gy = 1 degree) caused the vehicle 
to steady out at 137 degrees vice 159 degrees as shown in figure 34. This initial roll angle 
also caused a steady state roll angle of 17 degrees as shown in figure 35. In turn, this 


Steady state roll angle caused the steady state vaw velocity (r) shown in figure 36 (Le. 


CN 
tsi 


motion is no longer restricted to the vertical plane). Figure 37 shows a plot of z versus 
and indicates that the vehicle is taking a helical ascent as dicussed by Booth [Ref 2: 304 
305]. Therefore, the numerical integration solution resulting in a steady state pitch angle c 
159 degrees (figure 33) is unstable in the horizontal plane as predicted by the linearize 


dynamic response analysis. 
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Figure 33. Numerical Integration Solution for Angle of Pitch (6) when 


Center of Gravity is Aft Xcp = - 1%) 
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Figure 34. Numerical] Integration Solution for Angle of Pitch (6) when 


= - 1% and Initial Roll Angle (9,) is 1 degree 
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Figure 35. Numerical Integration Solution for Angle of Roll (9) when 


Xcop =- 1% and Initial Roll Angle (9) is 1 degree 
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Figure 37. Numerical Integration Solution for Heave Velocity (z) when 


Xcp = - 1% and Initial Roll Angle (g,) is 1 degree 
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Some question still remains with regards to thc measure of stabilty of the ‘inverted 
pendulum’ solutions predictcd by both the linearized dynamic response analysis and the 
numerical integration analysis. The linearized dynamic rcsponsc analysis predicts a stable 
solution for the case when center of gravity is placcd aft of center of buoyancy 
(Xop = - 1 %) and the dive planc anglc (0,) is - 7 degrces (figure 10). The corresponding 
steady statc value for pitch angle was 118 degrecs (figure 7). A random persistent roll 


disturbance (9,,) was added to the numerical integration model and the results arc shown in 


figure 38. The solid line indicates the results whcn a small disturbance is added (9, 


centcred about 0.1 degrees), the dashed linc indicates the results when a large disturbance 


is added (9, centered about 1.0 degrees). As expected, the large disturbance caused the 
vehicle to roll over as shown by the resulting angles of roll (9) in figure 39. Howcver, the 


vehicle continued to remain stable during smal] constant random disturbanccs in the 
inverted position as shown by the resulting angles of roll (g) in figure 40. This indicates 


that indeed these ‘Inverted pendulum’ solutions have a significant measure of stability. 
Pp 2 5 
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Figure 38. Numerical Integration Solution for Angle of Pitch (6) when 


che” 1% and a Persistent Roll Disturbance is Added 
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Figure 39. Numerical Integration Solution for Angle of Roll (9) when 


Xcp = - 1% and a Large Persistent Roll Disturbance is Added 
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Figure 40. 


VI. CONCLUSIONS AND RECOMMENDATIONS. 


The steady state analysis resulted in four possible solutions provided that vehicle 
motion was restricted to the vertical plane. Analyzing the dynamic stability using the 
steady state results as nominal points generally indicatcs that only one (if any) of the 
four possible solutions will be stable. There are a few cases where two solutions are 
stable, but these cases (certain non-zero bow plane deflection angles) appear to be the 


exception and not the rule. 


The dynamic stability characteristics of submersibles can be separated with respect to 
vertical plane motions (u,w,q,8) and horizontal plane motions (v,p,r,9). 


It is possible for submersibles to be dynamically stable with respect to vertical plane 
motion in the inverted (belly up) position during ascents (Inverted Pendulum’ 


Stabilization). 
‘Inverted Pendulum’ stabilization is also possible in the horizontal plane. 


Submersibles are able to maintain this invcrted orientation (i.e. ascend belly up 


without rolling over) even under some persistent roll excitation. 


As a recommendation, the dynamic stability analysis should be expanded to include 
the case where angle of roll is 180 degrees, and the casc whcre anglc of roll is non- 
zero (1.e. g neither equals zero nor 180 degrees). Analyzing the 9 =180 degrees cases 
will only involve changing a few signs with regards to trigonometric functions; 


however, analyzing the non-zero cases will require significant effort. 


Furthermore, identifying and characterizing different stability regions over a range of 


variations of the system parameters should be the matter of future research. 


APPENDIX A 


STIX DEGREE OF FREEDOM EQUATIONS OF MOTION 


Source: Smith, Crane, and Summcey [Reference 1:11-16] 


1. LONGITUDINAL (SURGE) EQUATION OF MOTION 


mu 4 - xg (q2 + 12) + vg (pq - 1) + 2G (pr + q) 
Xpp p2 + Xqq q2+X,, 12+ Apr pr 
+ Ay U+ Aw WG t+ Ayp Vp + Ave VF 
+ uq | Agds Os + Agdb db] + XrOp urOy 
+ | ee IV XO an db) 
+ u2(Xos05 ds? + Xobob Ob2 + Xordz Or) - (W - B) sin 8 


-) Ba 
+ U- prop 


2. LATERAL (SWAY) EQUATION OF MOTION 
m v+ur- wp + xq (pqtr)-vG (p2 +r2) + ZG (ar - p) = 
~Ypp + Yr + Ypq pat Yar oF 
Pay VF pup + iy Ut Vvg Yq + Iwo WP Yay 
i Neuve Vy + Vgn ue or 


“nose 4 9 (V+xr) .. 
- Cpy h(x} (v+xrj& +Cpz b(x} (w-xq}© » +———-" dx 
Xtal | UeflX) 


+ (W-B} cos 6 sing 


3. NORMAL (HEAVE) EQUATION OF MOTION 

m Ww -uq+ vp + X¢ (pr - q) + yc (qr + p) - zg (p? + q?) = 
Zaqq + Lop p- + L£pr pr + Lrr Ke 
+ Zw Ww + Zquq + LZyp vp + Zyr VI 
+ Ze uw + Zyy v2 + U2 (Zee O<+Z hb dp) 
Xnose - } (w-x 
~Cpy h(x) (v+xr}? +Cp,z b(x) (w-xq)* | net! dx 
; Uc glx) 


Atail 


- (W-B} cos 8 cos 


4. ROLL EQUATION OF MOTION 
Ix pt (Iz -ly} qr + Ixy (pr - q) - lyz (q2 : 12} 
- xz (pq +1) +m _yg (w - ug + vp)- 2G (v + ur - wp) = 
Kp p +ky r + Kpg pq + Kar qr 
ye ee Kp up + Ky ur + Kyg vq + Kwp wp + Kwr Wr 
+ Ky uv + Kyy vw, + yg W- yp B) cos 8 cos o 
- (zGW - zpB) cos 9 sind + ue Kprop 
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5. PITCH EQUATION OF MOTION 
ly q+ (Ix - ly] pr - Ixy (qr + p} + lyz (pq =) 
lee pe - r2)-m XG (Ww - ug + vp] - ZG (u - vr + wp) = 
‘Mg q + Mopp p- SF Mpr pr + Mrr Te 
+ Myw + Mg uq + Myp Vp + Myr VI 
+ | My) uw + Myvy v2 + u2 (Ms< Os + Mab dp} 
Xnose 4 on (W-X 
+ CDy h(x) (v+xr}* +Cpz bx) (w-xq}* | A x ax 
Atail | Me flX) 


-(xG W - xp B) cos 6 cos 9 - (7G W - zp B) sin 6 


6. YAW EQUATION OF MOTION 
ais (ly - Ix} Pq - Ixy Ip? ; q’) aly (pr + q) 
dele ar > pl =m XG (v + ur - wp} - yG (u-we+ wa) = 
-++- Nyy + Np up + Nr ur + Nvq Vq tNwp Wp + Nwr WT 
“+ Ny uv + Nyw VW + Nér u- op 
*nose | (v+xr) 
- CDy h(x) (v+xr}- + CDz b{x) (w-xq)? | ~ X dx 
Atail Vet) 


+(xG W- xp B}cos 6 sing + (yg W- ypB) sin 8 + u- Nop 


APPENDIX B 


ROTATION SEQUENCE AND EULER ANGLE RATES 


1. ROTATION SEQUENCE FOR 90, 8 AND ¥ 


Smith, Crane, and Summey [Reference 1:8] descibe the transition from body fixed to 


inertial reference frames as follows: 


Since the equations of motion are referred to an axis system that is fixed for the 
SDV (swimmer delivery vehicle), and thus translates and rotates with it, the 
Orientation and position of the moving body axis system relative to a fixed inertial 
reference system must be specified. The orientation of the body axis system with 


respect to the inertial reference system is defined by the standard Euler angles ¥ 
(yaw), 8 (pitch), and @ (roll). The rotation sequence from the inertial reference 


system to the body system is ¥ ,6 , and 9 as shown in Figure B1 taken from Smith, 
Crane, and Summey [Reference 1:18]. 


2. EULER ANGLE RATES FOR 0,0 AND ¥ 


The Euler angle rates used along with the six equations of motion (Appendix A) in 


order to completely determine the motion of the submersible were specified by Smith, 


Crane, and Summey [Reference 1:20] to be: 


p+qsin dtan@+rcos gtan 0 
Q = qcos 9 - rsing 
sin r COS © 


cos 8 cos 6 
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(1) Vehicle-Centered 
Gravity-Directec System 
parallel to inertial 
axis system. 





(3) Roll-Resolved Flight 
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paom «AY 2 py rota tion 
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angle &. 


(2) Horizontal Flight 
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bout Z through yaw angle w. 
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Figure BJ. Unit Sphere Development of Euler Angles 
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APPENDIX C 


STEADY STATE COMPUTER PROGRAM 


PROGRAM STEADY 
STEADY STATE SOLUTIONS IN THE VERTICAL PLANE 
DIVE PLANE VARIATION 


REAL L,MASS,1X,1Y,12Z,1XZ,1YZ,IXY,LAMBDA 

REAL KPDOT, KRDOT, KPQ,KOR, KVDOT, KP,KR,KVQ,KWP, KWR, KV 
REAL KVW,KPN, KDB 

REAL MQDOT,MPP,MPR,MRR,MWDOT,MO,MVP,MVR,MW,MVV,MDS 
REAL MDB,NDRB 

REAL NPDOT,NRDOT,NPQ,NQR,NVDOT,NP,NR,NVQ,NWP,NWR,NV 
REAL NVW,NDRS 


DIMENSTON X( 9) 8h 9 aoe Eee) 


GEOMETRIC PROPERTIES AND HYDRODYNAMIC COEFFICIENTS 


aQaAN 


PI =4.0*ATAN(1.0) 
WEIGHT=12000.0 
L = ted 25 
RHO = 1.94 
G = S258 
CDO = OS 0057 
MASS =WEIGHT/G 
ChZ =). 520. 5* RHO 
XWW = 1.710E-01*0.5*RHO*L**2 
XWDS = 4.600E-02*0.5*RHO*L**2 
XWDB = 9.660E-03*0.5*RHO*L**2 
XDSDS =-1.160E-02*0.5*RHO*L**2 
XDBDB =-8.070E-03*0.5*RHO*L**2 
CcDO = CD0*0.5*RHO*L**2 
ZW =-—-3.020E-01*0.5*RHO*L**2 
ZDS =-2.270E-02*0.5*RHO*L**2 
ZDB =-2.270E-02*0.5*RHO*L**2 
MW = 9,860E-02*0.5*RHO*L**3 
MDS =-1.113E-02*0.5*RHO*L**3 
MDB = 1.113E-02*0.5*RHO*L**3 
OPEN(21,NAME='’ST1.RES’,STATUS=' NEW’ ) 
OPEN(22,NAME='’ST2.RES’, STATUS='NEW’ ) 


OPEN(24,NAME=’ST4.RES’,STATUS=’NEW’ ) 


( 
( 
OPEN(23,NAME='ST3.RES’ , STATUS=’ NEW’ ) 
( 
OPEN(31,NAME='’COEF.DAT’ ,STATUS='NEWw’ ) 
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Cr Ore) 


Der ine The. LENGIE xX ,) BREADIE BR, AND HEIGHT HA TERMS 


Mel eee 10529 712 0 
Bde ve): (290 a8 7120 
a) =. SH 87.8712 -0 
rae) ped 66.3712 20 
sae Te 0 
X(6) = 83.2/12.0 
Ch ae 91.2/12.0 
X(8) = 96.7 712-0 
oy) 6] 10322712. 0 
Bey = 0.00/12.0 
BR.) os 8.24/12.0 
BR(3) = 19.76/12.0 
BEG) = 920286712 20 
BR(5) = 31.85/12.0 
BR(6) = 27.84/12.0 
BR(7) = 21.44/12.0 
BR(8) = 12.00/12.0 
BR(9) = 0.00/12.0 


COMPUTE AREA AND CENTROID 
CALL TRAP(9,BR,X,AREA) 

DO 9 I=1,9 
VEC1(1I)=X(1I)*BR(TI) 
CONTINUE 

CALL TRAP(9,VEC1,X,XAA) 
XA=XAA/AREA 


WRITE (*,1002) 

READ (%*,*) DSMIND,DSMAXD,IDS 
DSMIN=DSMIND*PI/180 
DSMAX=DSMAXD*PI/180 

WRITE (*,1001) 
READ (*,*) 
VRire { *, 1003) 
READ (%*,*) DELB 
DELB=DELB*WEIGHT/100.0 
WRITE (*,1004) 


RATIO 


READ (*,*) XGB 

XGB=XGB*L/100.0 

WRITE (*,1005) 

READ (*,*) ZGB 

WRITE (*,1006) 

READ (*,*) XB 

XB=XB*L/100.0 

WRITE (*,1007) 

READ (*,*) ZB 

Wepre (31. *): RATIO.) DELB.. XGB,.2G5, XB, ZB 


719 


DO=i l=) Gb s 


DS=DSMIN+(DSMAX-DSMIN)*(I-1)/(IDS-1) 
IF (DELB.EQ.0.0) DELB=0 7000001 

IF (ZGB.E0.0.0) ZCB POscguag 
DB=RATIO*DS 


PX =XGB*WEIGHT-XB*DELB 
PZ =ZGB*WEIGHT-ZB*DELB 
DEN =CDZ*AREA* ( PX+XA*DELB ) 


LAMBDA=MW* DELB-PX* ZW+PZ* (XWDS+RATIO*XWDB) *DS 

ALPHA =-PX*(ZDS+RATIO*ZDB) *DS—-PZ*CD0+PZ* (XDSDS+ 
RATIO*RATIO*XDBDB) *DS**2+DELB* (MDS+RATIO 
*MDB)*DS 

BETA =PZ*XWW 

LAMBDA=LAMBDA/DEN 

ALPHA =ALPHA /DEN 

BETA =BETA /DEN 


A = 1.0+BETA 
B = LAMBDA 
C = ALPHA 


DET= B*x*2-4.0*A*C 

IF (DETeLT.0.0) 1GO T1O— 

WP=(-B+SQRT(DET))/(2.0*A) 

YY=—-XWW*WP **2-(XWDS*DS+XWDB*DB) *WP 
—-(XDSDS*DS**2+XDBDB*DB**2)+CDO 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS(WP) 

IF (WP.LT.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS(WP) 

THETA=ATAN2(YY,XX) 

USQ=DELB*SIN( THETA) /YY 

THETA=THETA*180/PI 

DSD=DS*180/PI 

LF (USO. L020) SGeeioms 

IF (WP.GE.0.0) U= SQRT(USQ) 

IF (WP.LT.0.0) U=-SQRT(USQ) 

W=WP*U 

WRITE (21,*) DSD,THETA,U,W,WP 


WP=(-B-SQRT(DET))/(2.0*A) 

YY=—-XWW*WP**2-(XWDS*DS+XWDB*DB ) *WP 
—~(XDSDS*DS**2+XDBDB*DB**2)+CDO 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS (WP) 

IF (WP.LT.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS(WP) 

THETA=ATAN2 (YY,XX) 

USQ=DELB*SIN( THETA) /YY 
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DSD=DS*180/PI1 
THETA=THETA*180/PI1 

Poss OuiT 0.0) Go: fore 

IF (WP.LT.0.0) U=-SORT(USQ) 
IF (WP.GE.0.0) U= SQRT(USQ) 
WeWP*U 

WRITE (22,%*) DSD,THETA,U,W,WP 


2 A «= -1.0+BETA 

DET= B**2-4.0*A*C 

Pepe Tere 0)yGo Te 1 

WP=(-B+SQRT(DET))/(2.0*A) 

YY=—XWW*WP**2-(XWDS*DS+XWDB*DB) *WP 

& —(XDSDS*DS**2+XDBDB*DB**2)+CDO 

IF (WP.LT.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS (WP) 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS (WP) 

THETA=ATAN2 (YY, XX) 

USQ=DELB*SIN( THETA) /YY 

DSD=DS*180/PI 

THETA=THETA*1B80/PI 

PP (USO.nT. 0.0) GO TO 4 

Tew We.Ghs0.0). U=-SORT( USO) 

IF (WP.LT.0.0) U= SQRT(USQ) 

W=WP*U 

WRITE (23,*) DSD,THETA,U,W,WP 

4 WP=(-B-SORT(DET))/(2.0*A) 
YY=—XWW*WP**2-(XWDS*DS+XWDB*DB)*WP 
& ~(XDSDS*DS**2+XDBDB*DB**2)4+CD0 

IF (WP.LT.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS(WP) 

IF (WP.GE.0.0) XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS (WP) 

THETA=ATAN2 (YY,XX) 

USQ=DELB*SIN( THETA) /YY 

DSD=DS*180/PI 

THETA=THETA*180/PI 

ie fUuse).b1.0,0) Go To 1 

IF (WP.GE.0.0) U=-SORT(USQ) 

IF (WP.LT.0.0) U= SQRT(USQ) 

WeWP*U 

WRITE (24,*) DSD,THETA,U,W,WP 


is CONTINUE 


STOP 
POO FORMAT {' ENTER BOW PLANE TO DIVE PLANE RATIO’) 
1002 FORMAT /( END eR. tN, (162%. AND INCREMENTS iN 
& DS (degrees)’) 
1003 FORMAT ( ENTER DELB (%W)’) 
O04 FORWAT -(° ENTER.XGB (41) ) 


~~ ~= 


= ~ 


5 | 


1005 FORMAT (° “ENTER ZGE @Gfcer ae 
1006 PORMAT (° ENTER XB (2L) 
1007 FORMAT (’ ENTER ZB (feet)’) 


END 
c 
SUBROUTINE @TRAP(N, A, eo, Oo 
C 
c NUMERICAL INTEGRATION ROUTINE USING 
C THE TRAPEZOIDAL RULE 
€ 
DIMENSION A(1),B(1) 
Ni=N-1 
OUT=0.0 


DO 1 I=)7Ni 
OUTI=0.5*(A(I)+A(I+1))*(B(I+1)-B(1)) 
OUT =OUT+O0UTI 
IP CONTINUE 
RETURN 
END 


APPENDIX D 


LINEARIZED DYNAMIC STABILITY COMPUTER PROGRAM 


@ PROGRAM LINEARIZED DYNAMIC STABILITY 
G 10 Z0 30 40 
G@2375078970123450/690123256769012345676901234567890123456 
@ 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION L,MASS,IA,1X,1Y,12Z 
DOUBLE PRECISION KPDOT,KRDOT,KPQ,KQR,KVDOT,KP,KR, 
& KVQ,KWP,KWR,KV, 
& KVW,KPN,KDB,MQDOT,MPP,MPR,MRR,MWDOT,MOQ,MVP,MVR, 
& MW,MVV,MDS,MDB, 
& NDRB,NPDOT,NRDOT,NPQ,NOR,NVDOT,NP,NR,NVOQ,NWP,NWR, 
& NV,NVW,NDRS 
C 
DIMENSION OAT(4 ,4))51(4,4),BETAl (4), ALFRI(4),ALF1I1( 4) 
DIMENSION BB1(4,4),8B2(4,4),2221(4,4 
DIMENSION ee eee 
DIMENSION WR1(4),WR2(4),W11(4),W12(4 
DIMENSION Pee ees 
G 
@ GEOMETRIC PROPERTIES 
C 
PI = 4.D0*DATAN(1.D0) 
WEIGHT= 12000.0 
Oped — sO 
IDSs =e oO 0 
LZ — alle a 101 Oe 0 
1g = Le aS 
RHO = 1.94 
G = Bee 
EDO = O20:0s7 
MASS = WEIGHT/G 
CDZ = 0.5*0.5*RHO 
EDO. = CD0*0.5*RHO*L**2 
Cc 
@ SURGE HYDRODYNAMIC COEFFICIENTS 
é 
XPP = 7.030E-03*0.5*RHO*L**4 
XQQ =-1.470E-02*0.5*RHO*L**4 
XRR = 4.010E-03*0.5*RHO*L**4 
XPR = 7.640E-04*0.5*RHO*L*+4 
XUDOT =-7.580E-03*0.5*RHO*L**3 
XWO =-] -YZ0E-01*0.5* RHO L== 3 


moo 


GOO) 


QYdOA0D 


XVP =-3,240E-03*0.57RHO7E=-2 
XVF = ly890E—02* 02 Sano ise 
XODS: = 2761005 —07 “Ul series 
XODB =-2 6005-02 *0F = “FHOtrL =: 
XRDR ==8. TGUE=04*0e5"FHO- be 
XVV = 5. 2905-02405" RUO@ wie 
XWW = 1.710E-01*0.5*RHO*L**2 
XVDR = 1 7308-—05* 0S Resi 
XWDS = 4.600E-02*0.5*RHO*L**2 
XWDB = 9.660E-03*0.5*RHO*L**2 
XDSDS =-1.160E-02*0.5*RHO*L**2 
XDBDB =-8.070E-03*0.5*RHO*L**2 
XDRDR =-1.010E-02%*0.5*RHO*L**2 
XRES = CDO70.5*RHO*Ls42 


SWAY HYDRODYNAMIC CGEFFICIENTs 


YPDOT = 1.270E-04*0.5*RHO*L**4 
YRDOT = 1.240E-03*0.5*RHO*L**4 
BOO, =y4 1255-03 * 05 *hoO eis 
YOR =-6.510E-03*0.5*RHO*L**4 
YVDOT =-5.550E-02%*0.5*RHO*L**3 
me = 3. 0552-0570 55% RHO. iL +s > 
iE = 2.970E-02%*0.5*RHO*L**3 
LVAD = 2.360E-02%*0.5*RHO*L**3 
YWP = 2.3505B-01*0.5*REO*L**3 
YWR =-1.880E-02*0.5*RHO*L**3 
YV =-9.310E-02*0.5*RHO*L**2 
YVW = 6.840E-02*0. 5*REO*L< <2 
YDRS =42.270E=02*0. 5*RHO* LZ 
YDRB =#+2.270E-02*0.5*RHO*L**2 


HEAVE HYDRODYNANTC COErFICie as 


ZQDOT =-6.810E-03%0.5*RHO*L**4 
ZEPP = 1.270E-04%*0.5*RHO*L**4 
(pens: = 6.6/0E-0370 5*REC* ra 
ZRR =-7.350E-03*0.5*RHO*L**4 
ZWDOT =-2.430E-01%*0.5*RHO*L**3 
ZQ =—1 .s505-01*025*RHO* ie 
ZVP =-4.810E=-02*0.5*RHO*E*es 
ZVR = 4.550E-02%0 5 =kHO=0 3s 
ZW =-3.020E-01*0. 5*RHO*L**2Z 
2VV =-6.840E-02*0.5*RHO*L**2 
Z2DS =-2.270E—02*0 .5*RHO*L**2 
Z2DB =-2.2/05=0270. > He fe 


ROLL HYDRODYNANTC3 COE EE Lere ie 


KPDOT =-1.010E-03*0. 
KRDOT =-3.370E-05-0" 
KPQ =—6.930E-0s70- 


5* RHO*L* +s 
5*RHO* Ls 
Spats (On Iie eS 
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qdaAY 


mMYaANY 


KOR = ] 
KVDOT = 1 
KE =-] 
KR =—§8 
ARVO =—5 
KWP =-] 
KWR = 
KV = 3 
KVW =-] 
FITCH 

MOQDOT =-1 
Fee = 5 
he a=. 
MRR =-2 
MWDOT =-6. 
MO =-6 
MVP = l 
MVR = ] 
MW = 9 
MVV =-2 
NDS =—-]1 
MDB = 1 


,OS0E—02"0. 
-270EB-04*0. 
.100E-02*0. 
-410E-04*0. 
-LI1S5E-03*0. 
-2/0E-04%*0Q. 
-390E-02*0. 
0555-05" 0. 
~870E-01*0. 


.680E-02*0. 
.260E-05*0. 
.040E-03%*0, 
.860E-03*0. 


SLOE-02*0. 


.860E-02*0. 
~HEUB=03% 0%. 
. /30E-02*0. 
.860E-02*0. 
Po LUE U2 Un 
pa SE 02% 07: 
FiiSe=024%0. 


5*RHO*L**5 
5*RHO*L**4 
5*RHO*L** 4 
5*RHO*L** 4 
5*RHO*L** 4 
5*RHO*L** 4 
5*RHO*L** 4 
5*RHO*L* * 3 
S* RHO-Ls= 2 


BY DRODYNAMN TG (CORT EICIENT 


5*RHO*L**5 
5*RHO*L**5 
5*RHO*L**5 
5*RHO*L**D 
5*RHO*L** 4 
5*RHO*L** 4 
5*RHO*L** 4 
5*RHO*L**4 
5*RHO*L** 3 
5*RHO*L* * 3 
5*RHO*L** 3 
pRB O* Liers 


Povo no DRODY NAMIC “COEFFICIENTS 


IaaANYN 


NPDOT =-3.370E-05*0.5*RHO*L**5 
NRDOT =-3.400E-03*0.5*RHO*L**5 
Neo S=2, 1 0E—-02* 0, 5*RHO*L** 5 
NOR = 22. /5Ve-0s70, o*hAaOre.* & 
NVDOT = 1.240E-03*%*0.5*RHO*L**4 
Oe =-8.405E-04*0.5*RHO*L**4 
NR =-1.640E-02*0.5*RHO*L**4 
NVQ = eee ls Ol run) =i 4 
NWP =-1.750E-02*0.5*RHO*L**4 
NWR = 7.350E-03*0.5*RHO*L**4 
NV =-7.420E-03*0.5*RHO*L**3 
NVW =-2.670E-02*0.5*RHO*L**3 
WS: “=—l cP 3b—02*0.. s*RHOrL** 3 


NDRB =+1.113E-02*0.5*RHO*L**3 


BEeANE THE LENGIS > AND@BREADIB BR TERMS 


eer = Oo od 720 
X(2) = Se ee) era 6 
X(3) = HO. 371220 
X(4) = ~O 00374 2.0 
Kio}. = eae ae ee 
X(6) = Soe 7) 0 
a = Sipe alec. U 
X(8) = Oe ae 
X(9) = POs eZlZ a0 


YAN 


OG: 


BEY 1) = 0), CHCA, 13 
BR(2) = 8224/1220 
Bees) = ToS 6y7t 28 
BR(4) = 2osby 28 
Bri.) = 3128571 2..0 
BR(6) = 27 20474 a 
BR(7) = 2144 ee 
BR(8) = 123 VIG Pore 
BR(9) = O7.00 72 2720 


COMPUTE AREA, CENTROID, ANE MOMENT OF INERG 


CALL TRAP 
DO 9 I=l1 


(9, ,BR,xA,2REA) 
9 

=X(1I)*BR(I) 

=X( 1) 2VEC tee) 
as 

CALL TRAP(9,VECI x a 
XA=XAA/AREA 

CALL TRAP( 9; VEC2Z 54,12) 


WRITE (%,1001) 
READ (*,*) IRES 
OPEN(31,NAME='COEF.DAT’,STATUS='OLD’ ) 


READ(31,%) RATIO, DELS) GopeezZ Cera 


BUOY= WEIGHT + DELB 
XG=XB+XGB 
ZG=ZB+ZGB 


MASS MATRIX COEBPEICrEN Ts: 


Bl(1,1)= MASS - XUDOT 
B1(1,2)= 0.0 

B1l(1,3)= MASS*ZG 

Bilt 1,4)=- 070 

Bilt 2, 1) = 2020 

B1(2,2)= MASS - ZWDOT 
B1(2,3)=-(ZQDOT+MASS*XG) 
Bagc2.. 4) =— 06 

B1(3,1)= MASS*ZG 
B1(3,2)=-(MWDOT+MASS*XG) 
B1(3,3)= IY-MQDOT 

Ba (34 = nO ae 

B47 00 

Bie 4 2) = 00 

B1(4,3)= 0.0 

Bi(4,4)= 136 
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ZB 


MYA 


BZ =I X=KEDOT 

BZ el y=s0 20 

B2Zx ls) == (RV DOTSEMASS* ZG) 
B2(1,4)=-KRDOT 

eg) = OO 80 

7 aera ea 8. 
eo) 

Bae 4) = 0.0 

B2(3,1)==—( YPDOT+NASS*2ZG)} 
B2(3,72)= 0.0 

BZ 2,3) = MASS -YVDOT 
B2(3,4)= MASS*XG-YRDOT 
B2(4,1)=-NPDOT 

B24 72.)=— 020 

B2(4,3)= MASS*XG-NVDOT 
B2(4,4)= IZ-NRDOT 


OFEN(CES NAME="DEOS2Z. RES’ 7 STATUS=" NEW) 


Peo he See sul) GO” Om 
Ee Lae BO 02) GO SLO. AZ 
feeb Ss 2b Omo GO" LOS 


ie, CERES. EOe4 o GOTO <4 

OPENS (2 NAtVE= scllwREs ws talUs= “On” ) 
READ (21% ,END#=100) DSD,)THETO, U0,WU, WP 
GO TOS 

OPEN (22, NAME="ST25RES* ,STATUS=" OLD") 
READ 1227, 7) END=100).-DSD,;THETO, U0 TWO, WP 
GO TO 5 

OPEN (23,NAME= STS. RES -STATUS="OLD" } 
READ 23,7*, END=100) DSDPIREYO,U0,WO,WE 
GO TO. 5 

OPENV(UZ4,NANE="S14,RES ,STATUS=" OLD} 
Pee De 24,7, ENe-1o0) DSD,THETO U0, W0,WeP 
GO TO 5 


THETAO=THETO*PI/180.0 
DS = DSD*PI/180.0 
DB = DS * RATIO 


DAMPING NATRIXN .CORPPICIENTS 


Al(1,1)=-2.0*U0*CD0+W0?r (XWDS> DS+XWDB*DB ) 
+2.0*U0*(XDSDS”DS**2+XDBDB*DBE**2) 
1(1,2)= 2.0*XWW*W0+U0* (XWDS*DS+XWDB* DB) 
41(1,3)= (XWQ-MASS) *W0=-(XODS-DS+KXODB*DB)+U0 
1,4)=-(WEIGHT-BUOY)*DCOS(THETAQ) 


Oy) 


GZ 
fiat 


8 2 
Cul 


Al(2,1)= ZW*W0+2.0*U0* (ZDS*DS+ZDB*DB ) 
Al(2,2)= ZW*U0-2.0*CDZ*AREA*DABS (WO ) 

A1(2,3)= (ZQ+MASS)*U0+2.0*CDZ*AREA*XA*DABS (WO) 
A1(2,4)=-(WEIGHT-BUOY)*DSIN(THETAO) 

Al ( )= MW*wW0+2.0*U0* (MDS*DS+MDB*DB ) 


Saw 
Al(3,2)= MW*U0+2.0*CDZ*AREA*XA*DABS (WO ) 
3,3)= (MQ-MASS*XG) *U0-MASS*ZG*W0 
—-2.0*CDZ*IA*DABS (WO ) 
Al(3,4)= (XG*WEIGHT-XB* BUOY )*DSIN(THETAO)- 
(ZG*WEIGHT-ZB*BUOY ) *DCOS( THETAO ) 


Bi 4 = 020 

Ala e2 as 00 

Aa aS je 0 

Al(4,4)= 0.0 

A2(1,1)= KP*U0+(KWP-MASS*ZG) *W0 

A2(1,2)=-(ZG*WEIGHT-ZB*BUOY) *DCOS( THETAQ ) 

A2(1,3)= KV*U0+KVW*W0 

A2(1,4)= (KR+MASS*ZG) *U0+KWR*WO 

AZ(2,1)= 1.0 

DEA aA NM (0) 10) 

A2(2,3)= 0.0 

A2(2,4)= DTAN(THETAO) 

A2(3,1)= YP*U0+(YWP+MASS ) *w0 

A2(3,2)= (WEIGHT—BUGY )*DGOsS( fHETAY 

A2(3,3)= YV*U0+YVW*WO0-CDZ*AREA*DABS (WO ) 

A2(3,4)= YWR*W0+(YR-MASS ) *U0-CDZ*AREA* XA 
*DABS (WO) 

A2(4,1)= MASS*XG*wW0+NP*U0+NWP*WO 

A2(4,2)= (XG*WEIGHT-XB*BUOY) *DCOS( THETAQ ) 

A2(4,3)= NV*U0+NVW*WO0-CDZ*AREA*XA*DABS (WO ) 

A2(4,4)= (NR-MASS*XG) *U0+NWR*W0-CDZ*IA*DABS (WO ) 


BB Gy) Sayed) 
CONTINUE 
CONTINUE 


BE2t 2, 
CONTINGE 
CONTINUE 
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Ghiin RGG (4947s (Bol ALLER ALFIily SETA O. 2221718) 
CALL DEGSTB(DEOS1,ALFR1,ALFI1,BETA1,FREQ],WR1,WI1) 


GU rGet4 (at Serer Row ADP I? BETA 0, Zoe eT eR 
GALL DEGSTE( DEOS2 ALPR2.ALEI2 BETA? PREO? WRo- Wi?) 


IP eCDEOSI. GESDEQSZ) DEOSEDEOS! 
TP (DEOsd .LT DEOSZ |» DEOS=SDEOSZ 


WRITE (41,2001) DSD,THETO U0, K0, WR, DEOS 
& DEOS1 ,DEOS2 . 
Tes DEOS .LT0.D0 ) 
& WRITE (42752001) DS; THETO, UO, WC) WE. DEOS. 
‘ DEOS1 ,DEOS2 
Lh eV DEOsSi. UTC. DO) 
s WRITE (43,2001) DSD,THETO,U0,W0O,WP,DEOS, 
& DEOS1,DEOS2 


iSO - TO A 
ae Ce ec S ape 
GOTO? 
eGo Tor 14 


iF -(1TRES.EO. 2 
iP “{2PRES.EO% zZ 
iP (TREES. BO.Ss 
Pr {IT RES. EO.4 
G 
OC 25 TOP 
1001 FORMAT (’ ENTER THE RESPONSE DATA FILE DESIRED 
& awe ess SOR wee) wea) 
Z00}, “FORMAT (8515.5) 
2002 FORMAT (F10.3) 
END 
C 
SUBROUTINE DEGSTB(DEOS ,ALFR,ALFI,BETA, OMEGA,WR,WI) 
t(VPeie tT DOUBLE PRECISION (A=H,O-2Z) 
DIEMNENSTON- ALER(4 |} AGPI(4),BETA(4),wRt( 4), wil 2) 
Biel iid 
WR(1)=ALFR(1)/BETA(I) 
WI(I)=ALFI(1I)/BETA(I) 
iECONTINUE 
DEOS=-1.0E+10 
D@-2 5 i— 14 
LeRoi eDEOownGoOsrO 2 
DEOS=WR(I1) 
IJ=I 
2. GONTINUE 
OMEGA=W1(1IJ) 
OMEGA=DABS (OMEGA) 
RE TORN 
END 


SUBROUTINE TRAP(N,A,5,0UT) 
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AAANAN 


NUMERICAL INTEGRATION ROUTINE USING 
THE TRAPEZOIDAL RULE 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A(1),B(1i) 

Ni=N-1 

OUT=0.0 

1B) @ es NE I 
OUTI=0.5*(A(I)+A(I+1))*(B(I+1)-B(I)) 
OUT =OUT+O0UTI 

CONTINUE 

RETURN 

END 
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